function y = func_delta_2(q, m, CONSTS)

    k0 = CONSTS.k0;
    a = CONSTS.a;
    eta = CONSTS.eta;
    eps_a = CONSTS.eps_a;
    p = func_p_from_q(q, CONSTS);
    J_m1 = func_J_m1(q, m, CONSTS);
    J_m2 = func_J_m2(q, m, CONSTS);
    H_m2 = func_H_m2(q, m, CONSTS);
    J_m1_1 = func_J_m1_1(q, m, CONSTS);
    J_m2_1 = func_J_m2_1(q, m, CONSTS);
    n1 = func_n1_from_p(p, CONSTS);
    n2 = func_n2_from_p(p, CONSTS);

    y = n2.*((eta/eps_a).*J_m1.*J_m2 - (J_m1+(eta/eps_a)*J_m2).*H_m2) - ...
              n1.*((eta/eps_a).*J_m1_1.*J_m2_1 - (J_m2_1+(eta/eps_a)*J_m1_1).*H_m2) + ...
              (n2-n1).*(H_m2.^2-((p.^2)./eps_a).*((m^2)./((k0*a*q).^4))) + ...
              p.*(eta/eps_a).*(m./((k0*a*q).^2)).*(J_m1+J_m1_1-J_m2-J_m2_1);

    if(false)
        
        eps = CONSTS.eps;
        g = CONSTS.g;
        q1 = func_q1_from_p(p, CONSTS);
        q2 = func_q2_from_p(p, CONSTS);
        Q1 = k0*a*q1;
        Q2 = k0*a*q2;
        bessel_m_1 = besselj(m, Q1);
        bessel_mp1_1 = besselj(m+1, Q1);
        bessel_m_2 = besselj(m, Q2);
        bessel_mp1_2 = besselj(m+1, Q2);
        
        y_s = -1i*((eta-eps)/g)*(eta/eps_a)*(1./(((k0*a)^2)*(sqrt(eta/eps))*q)).*...
            ((bessel_mp1_1./bessel_m_1)-1i).*((bessel_mp1_2./bessel_m_2)-...
            (eps_a/eta)*(sqrt(eta/eps))*1i);
        
        figure; plot(q, real(y), 'b.-', q, real(y_s), 'r-'); title('Re(\Delta^{(2)})')
        figure; plot(q, imag(y), 'b.-', q, imag(y_s), 'r-'); title('Im(\Delta^{(2)})');
        
    end
    
end
    
    

    
    
    
